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The fully two-dimensional Peierls barrier map of screw dislocations in body-centered cubic (bcc) iron has been 
calculated using the first principles method to identify the migration path of a dislocation core. An efficient 
method to correct the effect of the finite size cell used in the first-principles method on the energy of a lattice 
defect was devised to determine the accurate barrier profile. We find that the migration path is close to a 
straight line that is confined in a {110} plane and the Peierls barrier profile is single humped. This result clarifies 
why the existing empirical potentials of bcc iron fail to predict the correct mobility path. A line tension model 
incorporating these first-principles calculation results is used to predict the kink activation energy to be 0.73 eV 
in agreement with experiment. 
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1. Introduction 

Recent nanopillar and nanowire experiments 
[T], as well as in situ observation of dislocation 
motions [2] , clearly demonstrate the highly pecu- 
liar natures of plasticity in body-centered cubic 
(bcc) metals. Plasticity in bcc metals is mainly 
mediated by the thermal activation of kink pairs 
in screw dislocation lines owing to strong lat- 
tice friction, and shows a strong dependence on 
the metal elements, the direction of the applied 
stress and temperature |3l4j . In face-centered cu- 
bic metals, the lattice friction is extremely weak 
and a dislocation motion can be described well 
by the universal model of phonon drag. There- 
fore, a reliable model of dislocation mobility in 
bcc metals is indispensable for simulating plastic 
deformation in such metals. 

One of the remarkable properties of a screw 
dislocation in bcc metals is that its motion is not 
confined in a plane but varies with the orienta- 
tion of the applied stress. An improper model 
leads not only to qualitatively incorrect mobility 
but also to a quantitatively wrong direction of 
migration. Fig. [T] shows several possible migra- 
tion paths for a screw dislocation which moves 



from one stable "easy-core" position to another. 
There are several unstable or metastable dislo- 
cation positions, which are referred to as "hard- 
core" and "split-core" positions. Middle points 
between easy-core positions (M points) are also 
shown in Fig. [TJ 

Anisotropic linear elasticity solutions (LES) for 
each dislocation configuration seen from the [111] 
direction are shown in Fig. [3] using the differen- 
tial displacement map (see Appendix A) . The 
positive and negative numbers in the figure rep- 
resents the helicity of the dislocation core. The 
solution for the M point is obtained by a simple 
interpolation between the solutions of two adja- 
cent easy-core configurations. This solution can 
be regarded as the sum of displacement fields in- 
duced by two "half" dislocation helicities located 
at the two easy-core positions. The solution for 
the split-core configuration can then be obtained 
by linear extrapolation using the hard-core con- 
figuration and the M point configuration. This 
solution can be regarded as the sum of displace- 
ment fields induced by two positive dislocations 
at the easy-core positions and one negative dis- 
location at the hard-core position. Note that the 
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LES for the easy-core and hard-core configura- 
tions have -D3 point group symmetry in Fig. [21 
as is also true for the results of density functional 
theory (DFT) calculations for bcc iron and other 
bcc metals ^SiGiT.Sj . We will later see that the dis- 
placement of the elasticity solutions for the split- 
core and M point configurations is also close to 
that of the DFT calculations. 

If either the hard-core or split-core configura- 
tion is metastable, the migration of a dislocation 
follows a bent line and the Peierls energy becomes 
double humped. On the other hand, when neither 
of the intermediate configurations is metastable, 
the path is close to a straight line and the Peierls 
energy becomes single humped. The difference in 
this energy landscape strongly affects the kink nu- 
cleation energy of a screw dislocation in different 
slip directions, and this determines not only the 
average migration velocity but also the average 
migration direction. Indeed, strong temperature 
and stress dependence on the average migration 
direction of a screw dislocation is reported in a 
molecular dynamics study of bcc iron screw dis- 
location [9]. 

However, since no quantum mechanical data 
of dislocation movement are used for the fitting 
of the empirical potential of bcc iron employed 
in Ref. j9j, the dislocation motion obtained in 
the molecular dynamics (MD) simulations should 
critically be checked by the DFT-based simula- 
tions. The Peierls energy in bcc iron without 
shear stress has been calculated by Ventelon and 
Willaime [10] using the DFT method with a lo- 
calized basis. In the present work we employ the 
more accurate plane-wave-based DFT to calcu- 
late a Peierls energy of an isolated screw dislo- 
cation core with and without the applied shear 
stress, obtain a complete Peierls barrier map for 
the two dimensional motion of a screw dislocation 
and identify its migration path. By coupling the 
DFT result with the line tension model of a dislo- 
cation, we also estimate the average slip direction 
of screw dislocations for various temperatures and 
applied stresses. 

This paper is organized as follows. We first de- 
scribe the method used for the DFT calculations. 
Further details of the DFT calculations to esti- 
mate the Peierls energy are dealt with in Section 



3. In Section 4 we describe a method to estimate 
and correct the finite size effect of DFT calcu- 
lations using the Green's function method. In 
Section 5 the results of the DFT calculations are 
summarized. Section 6 presents the line tension 
model of a dislocation to estimate the kink pair 
nucleation energy for various applied stresses. 
Concluding remarks are given in Section 7. 




Figure 1. Several possible migration paths of 
a screw dislocation in bcc metals. Black, gray, 
and white circles represent the dislocation cen- 
ter positions of easy-core, hard-core, and split- 
core configurations, respectively. Black squares 
arc middle points between easy-core positions (M 
points). Metal atoms are located at the white 
circles. 



2. The method of DFT calculations 

The electronic structure calculations and the 
structure relaxations by force minimizations in 
the DFT steps are performed using the Vienna ab 
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Figure 2. Differential displacement maps of the 
linear elasticity solution for easy-core, hard-core, 
M point, and split-core configurations. Filled cir- 
cles indicate the dislocation core position in each 
configuration; the positive and negative numbers 
indicate the dislocation helicity located at each 
triangle. See the main text for details. 



initio simulation package, with the projector aug- 
mented wave method and ultrasoft pseudopoten- 
tials. The exchange correlation energy is calcu- 
lated by the generalized gradient approximation. 
Spin-polarized calculations are employed in all 
cases. The Methfessel-Paxton smearing method 
with 0.1-eV width is used. Structural relaxation 
is terminated when the maximum force acting 
on the movable degrees of freedom becomes less 
than 10 meV/A. The error of energy that comes 
from this termination is estimated to be 2 meV 
by comparing the energy with the calculation for 
the smaller threshold of 3 meV/A. The cutoff en- 
ergy for the plane wave basis set is 400 eV, and the 
convergence of Peierls barrier energy with respect 
to the increasing cutoff is confirmed. The cutoff 
energy for the augmentation charges is 511.368 
eV. 



Periodic dislocation quadrupole configurations 
in a parallelepiped cell shown in Fig. [3] defined 
by the following three edges. 



ei = 

62 = 

63 = 



^^112 + ^lyViio + 1^111, 
Vlll 



(1) 



are used in all calculations, where v^i2 = 
ao[112]/3, viio - ao[110], Vm = ao[lll]/2, 
and flo — 2.833A is a lattice parameter esti- 
mated by a calculation of single BCC cell of 
[aoOO] X [OaoO] x [OOao] with 20x20x20 k-points. 
DFT calculations are mainly carried out on the 
Lx X Ly = 15 X 9 lattice , and the smaller and 
larger lattice sizes oi x Ly = 9x5, 9x7, 
15 X 7, 21x11 are also employed to estimate the fi- 
nite size effect. Hereafter we introduce the Carte- 
sian coordinates X, Y, and Z which are parallel to 
the [112], [110], and [111] directions, respectively. 
Two dislocations with opposite helicity are placed 
in the cell to form a periodic quadrupolar array. 

The reciprocal lattice vectors of the par- 
allelepiped described above are {k^, —ky/2,0), 
(0,fcj^,0), and {0,—ky/2,kz) where k^ = 
%/67r/(2aoLx), ky = \/27r/(2aoiy), and fc^ = 
2y37r/(3ao). Since the Monkhorst Pack k-point 
mesh for this inclined Brillouin zone breaks the 
symmetry of the independent inversions in the X, 
F, and Z directions, we use a Cartesian k-mesh 
of size kx/nx x ky/2ny x kz/lQ, which preserves 
the symmetry. The mesh number and Uy are 
adjusted so that the mesh size becomes close to 
0.lA~^. Convergence of the Peierls barrier en- 
ergy with respect to the increasing mesh number 
is confirmed. 

3. Calculation of the Peierls Barrier 

We denote the displacement of each atom rel- 
ative to the perfect crystal position by Ui — 
{u^,uf,uf) where i is an index of an atom. We 
also define a differentiated displacement as the 
normalized difference of the Z-component dis- 
placement for each pair of neighbor atoms Uij = 
{ul — Uj)/b, where b denotes the length of the 
Burgers vector. For both u and Uy , the min- 
imum image convention in the z axis is applied. 



and 



< 1/2 is always satisfied. We define 
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Figure 3. Periodic quadrupolar array of dislocation cores in a parallelepiped cell used in the DFT 
calculations, seen from the [111] direction. Gray circles are iron atoms and the two white circles are the 
two dislocations with opposite helicities. 



the two dimensional Peierls energy as the mini- 
mum energy under the constraint of dislocation 
location being at a certain position in the two di- 
mensional plane. By watching the differentiated 
displacement Uij , one can determine in which tri- 
angle a dislocation is located: a directed sum of 
three Uij around a triangle is equal to the dislo- 
cation helicity inside it. 

It is easy to show that a set of possible val- 
ues of three Uijs around a triangle which has 
non-zero dislocation helicity becomes an equilat- 
eral triangle in the three-dimensional space of 
Uij. Thus it would be natural to map this tri- 
angular region of u.^ to the triangle in the real 
space which surrounds the dislocation to define 
the intra-triangular dislocation position. Its ex- 
plicit definition is given as: 

P ^ -2{j1u32 + r^2Ui3 + flu2i), (2) 

where fl is the two-dimensional position of vertex 
i relative to the center of the triangle, as shown 
in Fig. m The position is relative to the center 
of the triangle. The Cartesian coordinate of the 
position is given as: 



For a given dislocation position {Px,Py)^ we 
define the Peierls energy Ep{Px, Py) as follows: 

Ep{P,,Py) = =minS(x^x°) 

= Eix^x^ix'^)) (3) 

where x'^ denotes a set of degree of freedom that 
is related to the dislocation position and is deter- 
mined by {Px,Py). The set of all the other de- 
grees of freedom is denoted by x°. In the present 
case, x''' are the Z displacements of the three 
atoms around a dislocation. The energy Ed{x'^) 
is a minimum of the total energy E{x'^,x°) for a 
given value of x"^, and x°{x'^) gives the values of 
x° which minimizes the energy. 

In the DFT calculations or MD simulations, 
Ed{x'^) can be easily obtained by fixing the values 
of and relaxing all the other degrees of free- 
dom. Moreover, one can calculate the derivative 
of Edix"^) with respect to one of the fixed degrees 



First-Principles Study on the Mobility of Screw Dislocations in BCC Iron 



5 



of freedom xf from the forces acting on xf which 
are also easily obtained, as follows: 



dEd^ 
dxf 



dE{x'^,x°{x'^)) 

dE{x'^, x") dE{x'^, x") dx°{x'^) 



dxf 



dx° 



dxf 



(4) 



where F^ is just a force acting on xf. The term 
dE{x'^, x°)/dx° is zero at x° because it gives min- 
imum of E{x'^,x°) with respect to x°. The force 
acting on the dislocation is calculated as follows: 



dEp 
dEp 



V2 



(2Ff - Fi ~ K 



3 )^ 



4{Fi 



Fi 



(5) 



It should be noted that the core position does 
not change when the three atoms move in the 
Z direction by the same amount d, and from 
the translational symmetry, Ed also remains un- 
changed. When there are two or more disloca- 
tions in the system separated by a long enough 
distance, energy difference induced by indepen- 
dent translations of Xd for each dislocation mainly 
comes from the long-range elastic interaction be- 
tween two dislocations, and is minimized by using 
relative relation of Xd, which is given by a linear 
elasticity solution. 

In general, Peierls energy Ep as a function 
of dislocation position P is discontinuous at the 
boundary between two triangles, since the choice 
of Xd differs in each triangle and different values 
of x° are obtained for different minimization con- 
dition. To deal with this discontinuity, we extend 
the mapping domain of Eq. in a hard-core 
triangle to include three adjacent easy-core posi- 
tions, as shown in Fig. [5] The values of Uij which 
correspond to these three extended vertexes be- 
comes Uij = ±1/3 on the boundary edge between 
the current hard-core and adjacent easy-core tri- 
angle and Uij = =f1/6 on the other two edges. 
(Note that when the dislocation position crosses 
over a boundary between two triangles, Uij on the 
boundary edge discontinuously changes from 1 /2 
to —1/2 if the minimum image convention is used. 
Instead, we allow Uij to become larger than 1/2 



atom 1 




atom 3 



atom 2 



Figure 4. Dislocation core position given by 
~2{fiU32 +'r2Ui3 + f3U2i) , whcrc fl is the position 
of vertex i relative to the center of the triangle. 
See the main text for details. 



outside the original triangle so that the position 
given by Eq. ([2]) becomes continuous.) These co- 
incide with the values of when the dislocation 
is at the center of the adjacent easy-core triangle 
if the easy-core state is symmetric as is the 
case of bcc iron. 

Therefore, by controlling the Z displacement 
of the three atoms around a hard-core triangle 
and relaxing all the other degrees of freedom, we 
can obtain a hard-core configuration, three adja- 
cent easy-core configurations, and any intermedi- 
ate configurations between them. A continuous 
energy landscape of an area which covers all the 
relevant configurations for a dislocation core mi- 
gration is also obtained. Later we will see that 
the split-core configuration has a very high en- 
ergy according to the first-principles calculation, 
in contrast to the MD empirical potential calcu- 
lation, and the migration path does not approach 
the split-core vertexes in the hexagon but always 
remains inside the hexagon. 

From the symmetry, the energy landscape in- 
side a hexagon can be determined by investigat- 
ing the energy inside a triangle the three vertexes 
of which are easy-core, hard-core, and split-core 
positions. Peierls energies for the intra-triangular 
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positions denoted by PO to P18 shown in Fig. [5] 
are calculated using the DFT, and the two dislo- 
cations are set to the same intra-triangular posi- 
tion. 




Split 



Figure 5. The domain of mapping between 
the core position and the Z-displacement of three 
atoms around a hard-core position is extended to 
a regular hexagon which includes three adjacent 
easy-core positions. 



4. Finite Size Effects 

Clouet et al. [TT] showed that the difference in 
the dislocation displacement field between LES 
and the DFT calculation can be described well 
by the line forces near the core, which come from 
the heavily displaced atoms at the core region. 
While this "core force" is very localized at the 
core region, the displacement owing to this force, 
referred to as the core field, is inversely propor- 
tional to the distance from the core, and the in- 
terference between the core fields induced by the 
multiple dislocations and their mirror images re- 
sults in finite size corrections of various quanti- 
ties. In the present work we directly observe the 



Hard 




Split 



Figure 6. The intra-triangular positions PO to 
P18 used in the DFT calculations. 



core force using the DFT calculation for configu- 
rations given by LES, and calculate the core field 
by the lattice Green's function method [S] to es- 
timate the finite size effect. 

Figure [7] shows forces acting on each atom, F[, 
calculated by DFT for the case of LES of the easy- 
core, hard-core, split-core and M position config- 
urations in the 21 x 11 system. The core force 
for the smaller case of 15 x 9 is very close to the 
21x11 case, and the maximum difference of forces 
between them is 0.02 eV/A. In the easy-core con- 
figuration, the distance between the innermost 
three atoms and the second innermost atoms is 
about 96% of that of a perfect crystal, and the 
second innermost atoms are pushed outward. In 
the hard-core configuration, the distance between 
the innermost three atoms is about 94% of that 
of a perfect crystal, and the innermost atoms are 
pushed outward more strongly than in the easy- 
core case. 

The DFT energy of the LES, -EfJJ, also shows 
a very weak size dependence. Let us consider the 
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relative value of E^^g with respect to the easy- 
core configuration divided by the number of dis- 
location, AE^^g . Hereafter the symbol A is also 
used for other quantities with the same meaning. 
This difference of energy mainly comes from the 
energy difference of core structure AEc and the 
difference of the interaction energy between dis- 
locations AEiNT- We assume that AE^^g = 
AEc + AEjpjT and Ej^jt is given by 



log(ri2) + Y Qiq j\og{ 

1=1,2 ] 



(6) 
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where /i = 57 GPa is the shear modulus of the 
strain in {110} < 111 > calculated by DFT, i and 
j denote the index of dislocation inside the cell 
and all the outside mirror images, respectively, 

= ±1 is the helicity of the dislocation i and 
is the distance between dislocation i and j. Fig. 
[5] shows dependence of AEc of the hard-core, M 
point and split-core configurations on the system 
size Lx y. Ly. Solid and dashed lines correspond 
to the data of the systems with the aspect ratio 
^/^Ly/Lx being greater and less than 1, respec- 
tively. These two cases show the opposite size 
dependence except for the split-core case in the 
9x7 system size. From this result, the finite size 
effect for AEc at the size 15 x 9 is estimated to 
be about 4 meV. 

The DFT energy of the relaxed configura- 
tion is denoted by E'^^^, and the energy dif- 
ference - £;f|J is denoted by 
This relaxation energy can be approximated by 
— i ^ ■ F^ul, where u1 denotes the change of dis- 
placement by the relaxation. This change of dis- 
placement can be approximately calculated by 
the linear elasticity theory, and in the present 
work we use the lattice Green's function method. 
We assume that a small perturbation 5 to a degree 
of freedom j induces a force HijS on the degree 
of freedom i, and that the Hessian matrix Hij, 
which is the inverse matrix of the lattice Green's 
function, is approximated well by that of a perfect 
crystal. Then the core field can be obtained 
by solving the following linear equations; 



Figure 8. The size dependence of the core energy 
AEc of the hard-core, M point and the split-core 
configurations. Solid and dashed lines correspond 
to the data of the systems with the aspect ratio 
\/?>Ly/Lx being greater or less than 1, respec- 
tively. 



To calculate Hij, a perfect crystal system is 
prepared using the same system size and bound- 
ary condition as those for the dislocation calcula- 
tion, and one atom is displaced in each of the X, 
F, and Z directions. The forces induced on each 
atom is then observed using the DFT calculation 
with the same parameters as the dislocation cal- 
culations. We find that the induced forces are 
proportional to the amount of displacement up to 
O.OSA, and that the force is only significant at the 
displaced atom and its nearest neighbors, while 
forces on the other atoms are negligible. If we de- 
note by H!^^ the matrix element for the displace- 
ment of the atom i in the direction a and the dis- 
placement of the atom j in the direction 6, the ob- 
served values are H^^^ = H^^ = -10.6 eV/A^, 
Hf^z ^ _Q_Q3 eV/12, and = for all a ^ b 
cases. For the neighbor atom j located at the +X 
direction of the atom i and displaced by —b/3 in 
the Z direction compared to the atom i, the ma- 
trix is i?^^ = 3.32 eV/A2, H^.^ = 0.201eV/l2, 



E 



ttZZ 



l.OOeV/1^ Hij^ = = 0.331eV/l2, 
and H,^^ = H^^ = Hf,^ = Hf^ = 0. Matrix 

''J '■J ''J ''J 

elements for other neighbor atoms can be calcu- 



(7) 
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lated from these values using the rotational sym- 
metry. The elastic constants expressed in the cu- 
bic axis derived from this matrix are Cn — 237 
GPa, Ci2 = 104 GPa, and C44 = 116 GPa, which 
are in reasonable agreement with the experimen- 
tally observed values Cn — 243, C12 = 145, and 
C44 = 116 GPa. The calculated anisotropy ratio 
2C44/ (Cii — C12) = 1.74 indicates that the elastic 
anisotropy is included in the matrix Hij , although 
the ratio is smaller than the experimental value 
2.36. 

The relaxation energy can then be estimated 
from the lattice Green's function as follows: 

-\h-^F^F';. (8) 



^RLX 



If E^[-^ agrees well with the actual value E^[^, 
one can estimate the relaxation energy of an iso- 
lated dislocation, E^[^ , using the core force of 
only one dislocation and using a sufficiently large 
lattice in the calculation of u^, and also estimate 



the finite size effect E^g^ 



^RLX 



T^GF 
^RLX ■ 



Finally, the Peierls energy lS.Ep for an isolated 
dislocation is estimated as lS.Ep — AE^^^ — 
AEipjx 



5. Results 

Tabic [1] summarizes energies of the various core 
positions shown in Fig. [51 calculated using the 
15 X 9 system. For the calculation of E^[j^, core 
forces smaller than 0.05 eV/A are omitted, and 
a periodic rectangle of 150 x 90 atoms system is 
used to calculate E^[^ . The agreement between 
ErEx ^"^^ Ej(lx is excellent, and we can use the 
Green's function method to estimate the finite 
size effect on the Peierls energy. Our estimate of 
the finite size effect is 2 meV for positions near the 
M point and the hard-core position, and 5 meV 
near the split-core position. Together with 
the 4 meV uncertainty of the finite size effect on 
AEc, the numerical error for the Peierls potential 
is estimated to be 6 meV for P9, P16, P17 and 
PIS, and 5 meV for the other cases. Fig. [9]shows 
a comparison of the core fields between the DFT 
calculations and Green's function method for the 
case of easy-core, hard-core, M point and split- 
core configurations. The agreement is excellent, 
except near the core region. Considering that the 



finite size effect mainly comes from the long-range 
part of the core field, we expect that the error 
of the core field near the core region calculated 
by the lattice Green's function method does not 
affect the calculation of finite size effect Epg^. 

The calculated values of the Peierls energy 
/S.Ep in Table [T] are fitted to a plane-wave ex- 
pansion which accounts for the periodicity and 
symmetry as follows: 

=EJC'l/e(^a) + C2/o(a;„)] 
+ Ea [C3/e(2x„) + C4/o(2x„)] 
+ [Ie[xi - X2) + !e{x2 " 2:3) + fe{xz - Xi) 

where Xa = f-ka, ka points toward three equiva- 
lent < 110 >directions with a length \ka \ = 
and Lq is the distance between adjacent easy-core 
positions, fe{x) = (1 — cos(a;))/2, and foix) = 
sin(a;)/2. The coefficients are Ci = 21.82, Ci = 
-14.51, Ca = 2.59, C4 = -2.72, and C5 = -2.89 
meV. Figure [TU] shows this fitted potential sur- 
face, with cross-symbols corresponding to the val- 
ues of AEp in Table [T] their relative values with 
respect to the fitted values are shown by perpen- 
dicular lines. 

The Peierls energy profile on a line segment be- 
tween the hard-core and the split-core position 
is shown in Fig. Ill[ which is a dividing ridge 
between two basins of easy-core potential min- 
ima. The minimum in this plot gives a Peierls 
barrier for dislocation migration, and one can see 
that the plot is nearly flat between the hard-core 
position and the M point. This indicates that 
the saddle point can be anywhere between these 
two points. However, energy minimization of a 
kink configuration will select the shortest path 
between the two easy-core positions and the ac- 
tual saddle point is expected to be near the M 
position. Fig. [11] also shows the energy pro- 
file calculated by two variants of the embedded 
atom method (EAM) potential by Mendelev et 
al. [13] , which differ from the DFT results signif- 
icantly. This difference results in widely different 
average slip directions, as will be demonstrated in 
later sections. The split-core position, which is 
metastable in the Mendelev potential, is actually 
unstable as has been shown in Ref. |12j . 

Our final estimate for the Peierls barrier is 35 ± 
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Table 1 

DFT energies and calculated correction terms for various core positions for the system size 15 x 9. All 
values are in the unit of meV. The numerical error for the Peierls energy AEp is estimated to be 6 meV 
for P9, P16, P17 and P18, and 5 meV for the other cases. See the main text for details. 
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5meV per 6, which is in good agreement with the 
experimental observation of 27 — 48 meV |2I14| . 
This result is also in good agreement with the 
DFT calculation of SIESTA code, 28 - 33 meV 

m- 

Next, we investigate the effect of external stress 
on the Peierls energy. A uniform shear strain eyz 
is applied to the system in DFT calculations by 
increasing the Z component of the cell vector e2 
by an amount ey^e^'i where is the Y com- 
ponent of e2 . This strain exerts forces in the X 
direction on the two dislocations. We investigate 
four cases, for the stain values of eyz = 0%, 0.5%, 
1.0% and 1.5%. When the two dislocations with 
opposite helicities move in the X direction by the 
same distance, change in the uniform strain com- 
ponent induced by this movement is canceled out 
and the uniform strain remains unchanged. We 



confirm that the variation in the uniform stress 
ayz is at most 2% in the DFT calculation when 
the two dislocations move from the easy-core po- 
sition to the M point. 

The effect of the applied stress on the Peierls 
energy is given by the Peach-Koehler term as fol- 
lows: 

{(cr* &) X f} • f, (10) 

where cr * is a tensor- vector product of the stress 
and Burgers vector, I is a vector parallel to the 
dislocation line and its length is equal to the dis- 
location segment under consideration, and r is a 
dislocation position. An additional effect of the 
applied stress, on the shape of the Peierls energy 
AEp, has been reported in Ref. [15]. It is im- 
portant to check the presence of this additional 
effect to model the behavior of dislocations under 
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Figure 10. Two-dimensional Peierls energy AEp, 
obtained by the DFT calculations. The curved 
surface is a plane-wave expansion fitted to the 
DFT results. The cross-symbols are the values of 
AEp in TabledJ their relative values with respect 
to the fitted values are shown by gray perpendic- 
ular lines. The thin vertical lines are numerical 
errors of DFT calculations. 



an applied stress. 

When the Burgers vectors of the two disloca- 
tions are anti-parallel, the term in Eq. (jlOp can- 
cels out and one cannot directly observe it. In- 
stead, we observe separately the forces and 
Fy acting on the two dislocations. By compar- 
ing the lattice stress Sx — —Fx/b'^± ayz and 
Sy = —Fy/lP' for the different strain cases, one 
can see whether the Peierls potential is affected 
by the applied stress or not. 

Fig. [12] shows a schematic of the method of 
calculation for the lattice stress. Since the ap- 
plied stress breaks the symmetry of inversion in 
the X direction, the lattice stress curve plotted 
against the reaction coordinate 77 can be asym- 
metric around ry — 1/2, as shown by the solid line 
in Fig. [121 For the "plus" dislocation, its position 
between the easy-core and the adjacent M point 
which is on the -\-X direction corresponds to the 
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Figure 11. Peierls energy profile of a line seg- 
ment between the hard-core and split-core posi- 
tions calculated by DFT. Dotted lines show cal- 
culation results obtained by two variants of the 
EAM potential in Ref.[13]. 



range < < 1/2, while for the "minus" disloca- 
tion it corresponds to 1/2 < 77 < 1. Observing the 
forces acting on the two dislocations separately, 
we can obtain the complete lattice stress curve. 

Figure [T^] shows a comparison of the lattice 
stress for different applied shear strains. Ob- 
served values of the uniform stress <tyz which 
are subtracted from the force for each strain case 
are also shown by horizontal lines. Although the 
curves do not conform to a single line, the devi- 
ation of Sx is best described by shifting to the 
—rj direction as the strain is increased, while the 
shape of the curves of Sy changes slightly for 
strains greater than 1.0%. If the applied stress 
only shift the curve and does not change its shape, 
the shape of AEp remains unchanged. Fig. [14] 
shows a comparison of the Peierls energy profile 
of a line segment between the hard-core and split- 
core positions for the eyz = and 1.5 % cases. 
The difference between the two cases is less than 
5 meV in the saddle point region. 

From these results, we conclude that the ef- 
fect of the applied stress is described well by 
the linear elasticity theory given by Eq. [TUj and 
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Lattice Stress 




Figure 12. Schematic of calculation of lattice 
stress curve from the forces acting on the two 
dislocations with opposite helicity. Black circles 
denote the dislocation positions where the lattice 
stress is observed in the DFT calculations. Solid 
and dotted curves represent the lattice stress of 
the zero stress and non-zero stress cases, respec- 
tively. Both dislocations are moved in the same 
direction, and the forces on dislocations with the 
positive and negative helicities are used to cal- 
culate the lattice stress for < ij < 1/2 and 
1/2 < ?7 < 1, respectively. 



its additional effect on the shape of Ep is very 
weak, at least up to 900 MPa, which is the max- 
imum stress investigated in the present work. 
The Peierls stress, above which the Peierls en- 
ergy landscape has no stable minima, is estimated 
from Eq. (0) to be 1000 ± 50 MPa. This stress 
is far stronger than the experimental estimation 
of ^ 400 MPa (such a discrepancy between atom- 
istic models and experiment is a general tendency 
in bcc metals |16|17ll8ll9j l. The Peierls stress is 
determined by the maximum slope of the poten- 
tial energy, and one can estimate its lower bound 
by AE / Axb~'^ , where AE is the saddle point en- 
ergy and Ax is the distance between the easy-core 
position and the saddle point. Since the value 
of AE agrees well with the experimental estima- 
tion and the lower bound Peierls stress calculated 




0.25 0.5 0.75 1 

Reaction Coordinate 



Figure 13. Comparison of the lattice stress for 
different applied shear strains of ey^ = 0%, 0.5%, 
1.0% and 1.5%. The uniform stress ayz observed 
in the DFT simulation for each strain is shown by 
horizontal dashed lines. 



with AE = 30 meV gives 690 MPa, which is 
still far stronger than the experimental estima- 
tion, this discrepancy of the Peierls stress cannot 
be attributed to the shape of the potential en- 
ergy landscape but should be attributed to other 
phenomena, such as interactions between disloca- 
tions |16) . 

6. Line Tension Model 

The dislocation mobility is determined by the 
formation enthalpy of a dislocation kink, a defect 
at which a dislocation line moves from one Peierls 
energy minimum to the next [17j . Fig. [15] depicts 
the kink nucleation mechanism of dislocation mi- 
gration. As shown in the inset of the figure, the 
total dislocation enthalpy reaches a maximum 
when some part of the dislocation line bulges and 
overcomes the Peierls barrier. After that, a fully 
formed kink pair moves in the opposite direction, 
lowering the total enthalpy as they separate un- 
til the whole dislocation line moves to the next 
minimum of the potential energy. 

The kink formation enthalpy calculated from 
the molecular statics [2D] or line-tension model 
adjusted on atomistic calculations |15I21| has 
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split-core position for eyz = and 1.5 % cases. 



been shown to reproduce the dislocation veloc- 
ity in the molecular dynamics simulations. For a 
screw dislocation in bcc metals, the kink width is 
estimated as 10 — 206 [22], and a direct calcula- 
tion of kink formation enthalpy requires at least 
several tens of slabs which amount to thousands 
of atoms. Since plane-wave basis DFT of thou- 
sands of atoms is not feasible, we use the following 
line tension model to estimate the kink formation 
enthalpy: 

Elt = -f J2j {Pj ^ Pj+iY 

+ Y.,^Ep{P,) + {{o*b)xl}-P,, (11) 

where K is a constant related to the elastic con- 
stants of the material, j is an index of thin 
slabs with the thickness h parallel to the Z axis, 
Pj — {Pj, Pf) is the position of a dislocation core 
in the slab j, AiJp is the Peierls barrier per Burg- 
ers vector h obtained in the previous section, and 
the third term is the contribution from the exter- 
nal stress. 

The constant K can be calculated by comparing 
the quadratic expansion of Eq. [11] in and 
with the expansion of Ed{x'^ + '5) hi (5, which is 



Figure 15. Schematic of dislocation movement 
by the kink nucleation mechanism, (a) A straight 
dislocation lies in a valley of Peierls potential en- 
ergy, (b) Some part of the dislocation bulges and 
overcomes the Peierls barrier, (c) A kink pair 
is fully formed and two kinks move in the op- 
posite direction, lowering the total enthalpy as 
they separate, (d) Both kinks are absorbed at 
the boundary of the dislocation, completing the 
movement process. 



given as follows: 

Ed{x'' + 5)- Edix^) ^ -J2S,F, + ^J2 H^A6,. (12) 

i ij 

The perturbation in x'^ also changes x°{x'^), but 
the change in x° does not affect the energy. The 
matrix Hij can be calculated from the change in 
the forces acting on each x'^ as follows: 

F,{x'' + 6)-F,{x'') = Y,H^J^J. (13) 

j 

First we demonstrate the calculation scheme of 
K using an EAM potential, and compare the kink 
shape predicted by the line tension model with 
the actual kink shape obtained in MD simulation. 
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An isolated screw dislocation of length 606 is pro- 
vided using the system consisting of 1452 atoms 
per layer. With a Mendelev potential of the vari- 
ant 5 [12] , this system is fully relaxed to the easy- 
core configuration while the outermost atoms are 
fixed to the linear elasticity solutions. Then one 
atom in the innermost three atom columns is dis- 
placed in the -\-Z direction by 0.01 A, and a force 
acting on every atom is calculated using the EAM 
potential. The induced force on the two atoms 
which are located directly above and below the 
displaced atom, is found to be 29 meV/A in the 
Z direction, while that on the displaced atom it- 
self is —82 meV/A in the same direction. Forces 
induced on the other atoms are negligible. Thus 
Hij is calculated to be Ha = 8.2 eV/A^ for the 
diagonal element and Hij = — 2.9 eV/A^ for the 
neighboring pair on the same atomic row. If we 
denote the Z displacement of the innermost three 
atoms on the layer j by Zj^i, Zj^2 and Zj^^, the total 
energy can be expanded in these quantities as 

j a=l,2,3 

with A = 2.4 eV/Aand Kq = 2.9 eV/A. Together 
with Eqs. © and (HH), we get K = ^Kq = 0.816 
eV/A^ for the Mendelev potential. The Peierls 
energy landscape of the same potential is calcu- 
lated using a single-layer, isolated dislocation con- 
figuration, and fitted to the form of Eq. [HI The 
coefficients are Ci = 40.4, C2 = 66.9, C3 = -0.1, 
C4 = 2.3 and C5 = 0.6 meV. 

Fig. [16] shows a comparison of kink shapes of 
the MD simulation and the line tension model 
with these parameters. In addition to the good 
agreement in the kink shape in the x direction, 
which was shown in previous studies using the 
line tension model |15] , the kink shape in the y di- 
rection also shows excellent agreement. The kink 
pair energy calculated by the MD relaxation is 
0.78 eV, while the line tension model gives 0.69 
eV. It has been shown that the kink shape and 
energy is asymmetric between the two types of 
kinks in the MD calculations [22]. This means 
that there should be higher-order terms of the 
gradient Pj — Pj-i in the energy, and the differ- 
ence of the energy between the MD and the line 



tension model can be attributed to this term. As 
a whole, these results clearly demonstrate the va- 
lidity of the calculation scheme of K presented 
above. By a simple scaling argument, one can 
show that the kink formation energy is propor- 
tional to ^JKEp, and fortunately the uncertainty 
in K only modestly affect the estimate of kink 
formation energy. 



-0.: 









^-i-i-i r-r-F^ 1 infi Tfinsinn Mndsl n___ 














MC 








































































10 20 30 40 50 
Layer 



Figure 16. Comparison of the kink shapes of the 
MD simulation and the line tension model. Px 
and Py are normalized by the distance between 
two easy-core positions. The inset shows the two- 
dimensional kink shape in the XY plane for both 
cases. 



To estimate the value K in the DFT calcula- 
tion, a single-layer hexagonal system consisting of 
108 atoms, which contains a single easy-core con- 
figuration in the center, is prepared and relaxed, 
while the outermost 33 atoms are fixed. The ob- 
tained configuration planes are then stacked on 
each other to form a three-layer system, and one 
of the innermost atoms is displaced in the Z di- 
rection by O.OlA. The forces induced by this dis- 
placement are calculated using the DFT. From 
the MD result and the DFT calculation of the 
perfect crystal, we found that forces on the sec- 
ond or farther neighboring atoms can be ignored, 
and that a three-layer system is suffice to calcu- 
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late the matrix Hij. The force induced on the 
two atoms located directly above and below the 
displaced atom is found to be 0.0308 eV/A; the 
forces on the atoms in the other atomic column 
are negligible. From this result, we get Kq = 3.08 
eV/A^, and consequently if = 0.866 eV/A^ 

The kink pair nucleation enthalpy Ek and its 
dependence on the applied stress is calculated by 
applying the transition state search method [33] 
to the model Fig. [T71 shows the kink shape 

predicted from the line tension model, when there 
is no applied stress. The kink pair energy and the 
kink width are estimated to be 0.73 eV and 106, 
respectively, which are in good agreement with 
previous studies [32] and an experimental esti- 
mate [24'. Since the kink width is proportional 
to by^ K/Ep and Ep is much higher than that of 
the Mendelev potential, the kink width is much 
shorter compared to the Mendelev potential case. 



(bO.8 

ns 
c 
■a 

o0.6 

o 

O 

§0.4 
o 

(0 

a> 

^0.2 




Width 



Pxf 

e 
-d}-- 
d 





■H, 







□ 

--D- - 
□ 



10 



20 



30 
Layer 



40 



50 



60 



Figure 17. Kink shape obtained from the line 
tension model with the parameters given by the 
DFT calculations. Px and Py are normalized 
by the distance between two easy-core positions. 
The kink width is estimated to be 106 . The in- 
set shows the two dimensional kink shape in the 
XY plane. 



Fig. [TH] shows the dependence of Ek on the 
applied stress for three maximum resolved shear 



stress directions, [111](110), [111](112) (twin- 
ning), and [111](2TT) (anti- twinning) , together 
with experimentally observed values from Ref . [24] 
and their linear fit. For the case of stress in the 
(110) plane, the direction of a dislocation migra- 
tion coincides with the stress direction, while for 
the case of twinning and anti- twinning, the high 
Peierls barrier around the split-core position pro- 
hibits a direct migration in the {211} planes and 
the macroscopic slip in the {211} planes is real- 
ized by random successive migrations in the two 
equivalent {110} planes, assuming that the screw 
dislocation is terminated by either a free surface 
or a periodic boundary. A trace of this random 
migration has been experimentally observed as a 
wavy pattern on a free surface of a thin film [2]. 
When the screw dislocation is part of a loop, the 
edge part energetically favors a straight line and 
random migration is expected to be suppressed. 
The overall tendencies shown in Fig. [THjsuch as 
the non-Schmit behavior that a slip in the twin- 
ning direction is easier than in the anti-twinning 
direction, agrees with the experiment in Ref. I24j , 
except for the already discussed discrepancy of 
the Peierls stress. 

The average slip direction at a finite tempera- 
ture is calculated using the Brownian dynamics 
simulation: 



Pj{t + At) = Pj - AtVjELT + V2k^^fAtrfj, (15) 

where V j denotes a derivative with respect to Pj , 
T is the temperature, ks is the Boltzmann con- 
stant and fjj is a two-component Gaussian ran- 
dom variable whose average and variance are 
and 1, respectively. It should be noted that the 
time scale in Eq. (fT5|) does not correspond to the 
physical one, and any momentum effect is not in- 
corporated. However, Eq. ([T5|) is able to capture 
the entropic effect on the migration frequency and 
the interaction between kinks. A time step of 
At = 1.1 X 10~^ is used, which ensures that the 
change in Pj at each step is at most 3% of the 
distance between two easy-core positions, and the 
iteration is repeated 10^ times for each case. 

Fig. dni shows the dependence of the average 
slip direction on the stress and temperature, for 
the case of shear stress in the [111] (110) direc- 
tion. The angle is relative to the (110) plane and 
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CI, the kink nucleation is so rare that we cannot 
obtain sufRcicnt data. Above C2, nucleation in 
the +60° direction becomes comparable to that in 
the 0° direction and the average slip direction be- 
comes positive. When the applied stress is greater 
than SOOMPa, the kink nucleation in this direc- 
tion becomes unstable and eventually becomes a 
kink pair in the 0° degree direction, and the curve 
C2 becomes a vertical line at SOOMPa. Above C3, 
kink nucleation in the 0, +60 and —60° direction 
all becomes comparable and the average slip di- 
rection decreases as the temperature is increased. 



Figure 18. Stress dependence of the kink 
nucleation enthalpy for three stress directions, 
[111](110), [111](112) (twinning), and [111](211) 
(anti- twinning) . Experimental data from Ref. 
[24] are also shown as symbols, together with 
straight line fits. 



is defined as positive for the twinning direction. 
Similar plots have been obtained based on the 
MD simulations |9l25j . and the main difference 
between the present result and the MD simula- 
tions is the magnitude of the deviation from the 
(lIO) plane. For all the cases studied here, it 
is less than 2 degrees, in contrast to the MD 
results, in which the deviation is as large as 20 
degrees. Experimental observation of dislocation 
motion [2] indicates that the slip plane for the 
case of maximum resolved shear stress in the 
< 111 > (110) direction is primarily (110), and 
we conclude that the deviation of the slip direc- 
tion in MD is an artifact of the EAM empirical 
potential. The three lines CI, C2 and C3 in Fig. 
[T9l which separate the regions of different slip 
behavior, are given by equations T ^ E^/20kB, 
T = {Et ~ El)/2kB and T = {E^_-~ El)/2kB. 
respectively, where E^, E'^ , and E/, denote the 
kink nucleation enthalpy to the neighbor easy- 
core configuration located in the directions of 0, 
+60 and —60°, respectively, for a given value of 
applied shear stress in the [lll](ll0) direction 
calculated using the line tension model. Below 
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Figure 19. Temperature and stress dependence 
of the average slip direction, relative to the maxi- 
mum resolved shear stress direction of [111] (Oil). 
The curves CI, C2 and C3 are the boundaries be- 
tween different slip behavior regions, calculated 
from the difference between kink nucleation ener- 
gies in different directions. See the main text for 
details. 



Fig. [201 shows how the applied stress affect 
the nucleation energies E'^ , and i^^T- In the 
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present study, migration paths are close to the 
straight hne and the contribution of the apphed 
stress a to the Peierls barrier is approximately 
proportional to acos{9)x, where < x < 1 de- 
notes a reaction coordinate of the dislocation mi- 
gration and 6 = 0, +60 and —60°. The migra- 
tion path for the +60° case can bend toward the 
energetically favorable direction and the migra- 
tion energy is further lowered, while for the —60° 
case the high energy of the split-core position pre- 
vents the path from bending to the energetically 
favorable direction and the migration energy is 
slightly higher than the +60° case. Thus we have 
< K < ^k' in which case the stress and 
temperature dependence of the slip direction in 
Fig. [19] is reproduced. 

For the case of empirical EAM potential, the 
migration path passes though the split-core po- 
sitions and the effect of the applied stress dif- 
fers significantly from the present study case, as 
shown in Fig. [20] (b). The migration energes 
for the 0° and —60° cases are almost the same, 
since the saddle point of the migration for these 
two cases are close together even though the cor- 
responding atomic configurations differ. Conse- 
quently, the average slip direction becomes close 
to —30° as soon as the temperature is high enough 
to overcome the small energy difference between 
the 0° and —60° migrations. The double-humped 
potential shape also causes discontinuity in the 
kink nucleation enthalpy at a certain applied 
stress [26] . 

7. Summary and Conclusion 

Essential properties of a screw dislocation in 
bcc iron have been calculated and determined 
by the plane-wave-based DFT calculations. The 
properties of dislocation kinks are determined us- 
ing a line tension model which is based on the 
DFT results. Our findings are summarized as fol- 
lows: (i) a screw dislocation migrates from one 
easy core position to an adjacent easy core posi- 
tion, following a nearly straight path; (ii) the sad- 
dle point of the migration is around a midpoint 
of two easy core positions, at which the Peierls 
barrier is estimated to be 30 — 40 meV; (iii) a 
kink pair nucleation enthalpy is estimated to be 
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Figure 20. Schematic of the effect of the applied 
stress on the average slip direction for the cases of 
(a) the present study and (b) the empirical EAM 
potential. Solid and dotted lines show migration 
pathways with and without the applied stress, re- 
spectively, (c and d) Schematic behavior of the 
Peierls barrier for both cases. See the main text 
for details. 



0.73 eV when no external stress is applied; (iv) 
when the maximum resolved shear stress is in the 
< 111 > {110} direction, the average slip direc- 
tion remains in the {110} plane with a deviation 
of less than 2 degrees; (v) when the maximum 
resolved shear stress is in the < 111 > {211} di- 
rection, two equivalent {110} slips are activated 
randomly; and (vi) the Peierls stress is estimated 
to be 1000 MPa, with a lower bound of 620 MPa. 

These results are consistent with experimental 
estimation, except for the far larger value of the 
Peierls stress (vi) compared to the experimental 
estimate. This discrepancy in the Peierls stress 
between atomistic scale calculations and experi- 
ments is a general tendency in bcc metals, and 
we foresee that some mesoscopic modeling is re- 
quired, such as the pile-up model proposed by 
Groger et al.[T2], to relate the calculated Peierls 
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stress to the experimentally estimated yield stress 
in the zero temperature limit. Points (iv) and (v) 
are in clear contrast to the molecular dynamics 
simulations of bcc iron [3]. This is because the 
fundamental properties of migration for a screw 
dislocation have not been included in the fitting 
targets in the previous development of empirical 
potential in bcc metals. Recently, an improved 
EAM potential has been developed by Gordon ct 
al. [57] in an effort to reproduce both the sym- 
metric easy-core structure and the single-humped 
Peierls energy. Gordon et al. concluded that 
conventional EAM potential is not capable of re- 
producing these two properties simultaneously. 
Further improvements in EAM potential based 
on the dislocation migration properties obtained 
in the present work, possibly by including the 
anisotropic terms or the magnetic degrees of free- 
dom [28l in the potential, are strongly expected 
for the reliable simulation of plastic deformation 
in bcc iron. 
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Appendix A: Anisotropic elasticity 

Owing to the elastic anisotropy, the elastic 
constants in the Cartesian coordinate in Fig. 
|3] change from the cubic axis case by amounts 

^Cyy^zz = ^Czz,xx = / 3 , lS.Cxx,yy — A/6, 
^Czx,ZX = ^Cyz.yz = A/3, lS.Cxy,xy — A/6, 

^Cxx,zx — ^Cxy.yz — 3y^' ^^yy^^^ 

where A = Cn — C12 — 2C44 is the anisotropy 
parameter. The isotropic linear elasticity solution 
for the < 111 > screw dislocation has component 
Cyz and Cxz, and induces excess stress field axx — 
-A'ezx, cTyy = A'szx, and axy = -A'eyz where 
A' = -A=. However, these excess stress fields sat- 
isfy the equilibrium condition J2b (^bO'ab — and 
the solution itself is identical to the isotropic case. 
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Figure 7. Core forces of the easy-core, hard-core, spht-core and M position configurations. Z-components 
of the force are shown by vertical arrows and exaggerated five times compared to the XY components. 
Forces less than 0.05 eV/A are omitted for clarity. 
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Figure 9. A comparison of the core fields calculated by the DFT (solid red arrows) and the lattice 
Green's function method (dashed blue arrows), respectively. Only the XY components are shown for 
clarity, and they are exaggerated by a factor of 40. 



